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Abstract. We introduce the quantum stochastic differential equation (QSDE) 
approach to exactly analyze the response of quantum systems to a continuous-mode 
two-photon input. The QSDE description of the two-photon process allows us to 
integrate the input-output analysis with the quantum network theory, and so the 
analytical computability of the output state of a general quantum system can be 
addressed within this framework. We show that the time-domain two-photon output 
states can be exactly calculated for a large class of quantum systems including passive 
linear networks, optomechanical oscillators and two-level emitter in waveguide systems. 
In particular, we propose to utilize the results for the exact simulation of the stimulated 
emission as well as the study of the scattering of two-mode photon wave packets. 


1. Introduction 

The study of the dynamics of photon-photon interaction is fundamental in physics. 
Particularly, the controllable photon-photon interaction may play a vital role in the 
realization of all-optical circuits and quantum information processing [I1121I3]- For 
example, the transmission of single-photon signal might be controlled by a gated photon, 
leading to a novel design of photonic transistor which operates with minimum energy 
usage pp. Since photons rarely interact in free space, quantum systems are often 
employed to mediate the interaction. There exist numerous proposals for the mediation 
of light-light interaction using quantum systems such as artificial atoms in waveguides 
and molecules [H El E]- 
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In the language of system theory, the mediated two-photon interaction can be 
understood as the response of a quantum control system to a two-photon input. The 
output state carries the full information of the response of the system, which can be 
further used for photon statistics and correlation analysis. The exact calculation of 
the two-photon response of a two-level system has been studied using input-output 
formalism [7] , Bethe-ansatz method |1] and Lehmann-Symanzik-Zimmermann reduction 
[H]. The generalized treatment of photon-photon interaction has also been studied in 
[HI Ho], focusing on the analytical property of the scattering matrix. A diagrammatic 
approach is studied in m which took into account the effect of the relaxation of 
two distant qubits in scattering. In general, these methods model the interaction 
between the photons and the system as an inelastic scattering process, which could 
facilitate the stationary state analysis in either frequency or time domain [T21 fTB] . 

The multi-photon response of linear (harmonic oscillators) and finite-level systems has 
also been investigated in m HSI using QSDE equations. These works have shown 
that some important quantities such as output photon flux and covariance function 
can be conveniently calculated using the QSDE equations. Generally speaking, it is 
straightforward to study the time-domain dynamics by solving the Schrodinger equation 
which governs the interaction between the photon wave packets and the control system 
[ISIIIZI. For example, numerical results for the scattering of photon wave packets by 
a two-level emitter in one-dimensional waveguide have been obtained using this wave 
function approach HU. Despite the progress, however, the exact calculation of the time- 
domain output state is still challenging due to the complexity of two-photon dynamics. 

In this paper we propose a QSDE approach for the general modelling of the 
two-photon process. QSDE [TB] is the generalization of the equation of motion in 
Heisenberg picture, which can characterize the evolution of the operators for a general 
quantum system or a coherent quantum network. Specifically, QSDE can be derived 
using the parameters of the overall system, while the Hamiltonian and field coupling 
operators of the overall system are calculated according to the interconnection between 
the subsystems HU HHl in]. Therefore, the QSDE approach is capable of dealing 
with an integrated quantum system which may involve arbitrary number of subsystems 

[231123] • Moreover, it is convenient to utilize the QSDE approach to calculate the 
analytical form of the time-domain output state. As we will show in Section]^ the time- 
domain output state can be analytically calculated for a general passive linear network. 

The exact two-photon response of an optomechanical system and a two-level emitter 
can be exactly calculated as well. 

The QSDE formalism for dealing with two-photon response of a general quantum 
system is developed in Section The analytical computability of the output state is 
discussed in Section]^ In Sectionsandwe show some applications of the QSDE 
approach. Particularly, we consider the simulation of the stimulated emission and the 
scattering of two photons in a waveguide, based on the exact calculation of the output 
state. Conclusion is presented in Section 
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Figure 1. A quantum system may have multiple input channels. The Hamiltonian 
and the coupling operators of the total system are determined by its internal structure, 
which may involve interconnection of subsystems. 


2. The input-output formalism for uncorrelated and general two-photon 
input states 

In this section, we propose a formalism to model the interaction between a continuous¬ 
mode two-photon input and a general quantum system. 

First, we introduce the QSDE description of the dynamics of open quantum systems. 
The dynamics of an open quantum system interacting with the input fields is generated 
by a unitary evolution characterized by a unitary operator U{t, to), where to is the initial 
time of the interaction. The dynamical equation of f/(t,to) is given by 

dU (t, to) = {b\t)L - L^b{t) - {^L^L + iHo)}U (t, to)dt, t > to, ( 1 ) 

with U{t + dt, to) = U (t, to) -b dU (t, to) and U (to, to) = / ® I- The system is coupled to 
the external fields through K channels (Figure]^, with the coupling operator defined 
by L = [Li ■ ■ ■ Lk]'^ which is a column vector of operators. Hq is the system 
Hamiltonian. b{t) = [6i(t) ■ ■ ■ bxit)]^ is a column vector of bosonic field annihilation 
operators, and Li,bi{t),i = 1,2, ■ ■ ■, K are defined on the t-th channel. The singular 
field operators satisfy the commutation relation [&i(t), &j(s)] = 6{t — s), i = j and 
[bi{t),b]{s)] =0, i^ j. Formally, Bi{t) = f* bi(s)ds is the quantum Wiener process and 
dBiit) = Biit + dt) — Biit) is the operator-valued Ito increment. Eq. ([^ is obtained 
by modelling the environment Hamiltonian as j'^^u:U{u:)b{u:)duj, and the interaction 
Hamiltonian as ff^(L'^b(u) -|- U(u)L)du; with L being independent of u. Here rotating 
wave and Markov approximations have been invoked to obtain the current solvable 
form of interaction Hamiltonian. It is worth mentioning that the rotating wave and 
Markov approximations are generally valid for quantum photonic systems mm, where 
interaction strength is relatively low, the incident photons are near resonance with the 
system transition frequency and the environment has no memory effects. The systems 
considered in this paper can be operated within this regime. 

Please also note that the choice of stochastic calculus would not affect the 
calculations of physical quantities [2S]. 
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Denote |0) as the vacuum field state and |0s) as the ground state of the system. 
Throughout this paper we concern with systems possessing the simple passivity property 

t/(t,to)|00,) = |00,), (2) 

by which we can prove 

f/^(t,to)|OOs) = f/^(t,to)f^(t,to)|OOs) = |00s). (3) 

The condition Eq. (|^ can be easily established as long as Hq and the couplings add 
no energy to the overall system. For later use, we make the following assumption 
throughout this paper: 

Assumption 1 The number of quanta in the overall system is a conserved quantity as 
time evolves. 


Under this assumption, the system will remain at the ground state when initially the 
system is at the ground state and the field is vacuum. The Heisenberg-picture evolution 
of a system operator X is defined by X(t) = W(t,to){I ® X)U(t,tQ), with / being the 
identity operator on the field. Based upon Eq. Q, the dynamical equations of X(t) are 
derived as QSDEs |26l |23 [HI l2H]: 

Xit) = C\Xit)) + b\t)[XitfLit)] + [Lt(t),X(t)]6(t), (4) 

bout{t) = L{t)+ b{t), (5) 

where the generator £l(X(t)) is given by 

£t(A'(i)) := -i[X{t),Ho(t)] 

« 1 1 
+ - Mt)nMx{t) - -xmaww))- (O 

k=l 

It is clear from Eq. Q that the evolution of a system operator is driven by the input 
field. Moreover, the output field operator boutit) is related to the input field operator 
b{t) by the following relation |2H] 

bout{t) = U\t,to)b{t)U{t,to) = U\T,to)Kt)U{T,to), T>t. (7) 


That is, the unitary evolution transforms b{t) into bout{t) in an infinitesimal time interval, 
which is a consequence of the Markov approximation. 

Next, we will define a continuous-mode two-photon state. Heuristically, we consider 
a singular input state &i(ti)&j(t2)|0) (See Figure ^ which contains two impulses of single¬ 
photon inputs at ti and t 2 , respectively. When the initial state of the system is |0s), the 
unitary evolution of the system-field state |tk(t)) is given by 


|vl/(t)) = f/(t,to)&l(ti)& 5 (t 2 )| 00 ,). 


( 8 ) 


Here ti,t 2 G is required in order to guarantee that U{t,to) covers the effective 

interaction process and so |T(t)) is the output state. Since the continuous-mode single¬ 
photon pulse can be modelled as the superposition of single-photon impulses [2H], the 
general form of an uncorrelated two-photon input state can be defined as 

/ oo poo 

dti / dt2b\ti)^i{ti)b\t2)^2it2)\0), 

OO J —oo 


( 9 ) 
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U(t, to) 



— 


1 


— 00 to ti t2 t 00 


Figure 2. The two-photon input is modelled as two separated S pulses. The first 
photon enters and interacts with the system at t = ti, then the second photon enters 
at t = < 2 - The unitary operator 17(t,to) governs the entire process. 


where each ^q(-) = {■), q = 1,2 is the collection of pulse functions over 

the K channels for each single photon. The pulse functions satisfy the normalization 
condition = 1. According to the generalized dehnition Eq. (9), 

each single-photon input could be superposed over the K channels. Since the input 
state Eq. 0 is dehned on (— cxd, -|-cxd), we must let f —)■ cxd and R —)■ —oo in order to 
obtain the correct output state. The joint (system plus field) output state can thus be 
calculated by 

|T(cx))) = f/(cx), -cx))|ljil ^2 0^) 

/•OO /‘OO 

dt2U{oo, -oo)b\ti)^i{ti)b\t2)^2it2)\00s)- 


dti 


( 10 ) 


Here U{oo, —oo) is well-dehned due to Eq. ([^. We need to convert Eq. (10) to a 
computable form. Note that the output state |\h(cx))) is a two-photon state with two 
excitations in the field, according to Assumption Here we have ignored the component 
^iMlois), r e ( —oo, oo) in the output state, where |ls) is the system state containing 
one excitation. This component vanishes in the steady-state limit by taking to —>■ —oo. 
More details on the steady-state limit can be found in the next section. Hence, the basis 
vectors of |T(cx))) are 

bl{Tl)b]{T 2 )\ 00 s),Ti,T 2 e (-CX),CX)). (11) 

As a result, the output field state is calculated by 


fieid{oo)) = {0s\ / dti / dt 2 U{oo,-oo)b\ti)^i{ti)b\t 2 )^ 2 {t 2 )\ 00 s). 


( 12 ) 


Inserting the identity, i.e. J^^dTidT2b\{Ti)b^j{T2)\00s){00s\bi{Ti)bj{T2) onto 

the two-photon subspace, we obtain 

^ poo poo 

\'^ fieldioo)) = ^ / dTidT2^[j{Ti,T2)bl{Ti)b]{r2)\0), (13) 

i,j=ld-oo J-oo 

with the coefficients {^ij{Ti,T 2 )} of the basis vectors defined by 


^'ij(,n,T2) = {00s\bj{T2)bi{Ti) / dti / dt2U{oo,-oo)b\ti)^i{ti)b\t2)f2it2)\00s). (14) 












Exact Analysis of the Response of Quantum Systems to Two Photons Using a QSDE Approach^ 


According to this expression, \^ij{Ti,T 2 )\'^dTidT 2 is the probability of emitting the two 
photons to the i-th. and j-th channels of the output held during [ri,ri + dri) and 
[r 2 , T 2 +dT 2 ), respectively. The output held state |T/jeid(c>o)) is analytically computable if 
T 2 )} are analytically computable. Inserting the identity t/(oo, —oo)W(oo, — 00 ) 
helps simplify the expression, which will enable the analytical computability of T 2 ) 
to be studied by 

/ oo poo 

dti / (if 2 f^(oo,- 00 ) 17 ^( 00 ,-oo) 6 j(r 2 )f/(oo,- 00 ) 

■OO J —00 

xU\oo, -oo)bi{Ti)U{oo, -oo)b\ti)ii{ti)b\t 2 )i 2 {t 2 )\^^s) 


— (00s I / dti / df2&j,ont (Ti)M(ti)^i(ii)6t(t2)^2(i2)|00,), 


(15) 


where bi^out{-),bj^out{-) are the i-th and j-th components of bout{-), respectively. Here 
we have used the property Eq. ([^ and the input-output relation Eq. ([^. Note that 
the pulse function ^'ij{Ti,T 2 ) is symmetric with respect to ri = T 2 , which is due to the 
indistinguishability of photons. Making use of Eq. (|^, Eq. (15) can be decomposed into 
four terms 


4(H,T-2) = 


dti 


dt2[Oi + Ojj + Ojjj + Ojv], 


(16) 


with 


dj = {00s\bjiT2)bi{Ti)b^{ti)fi{ti)b^{t2)^2{t2)\00s), 

Oil = (00s|Tj(r2)6i(ri)6'^(fi)^i(fi)&'^(f2)6(f2)|00s), 

OlII = {00s\bj{T2)Li{Ti)b'^{ti)fi{ti)b\t2)^2{t2)\00s), 
div = (00s|Tj(r2)Li(ri)6^(fi)^i(ti)&^(f2)6(^2)|00s). (17) 

Thus far, we have developed the formalism for the calculation of the response of a 
general quantum system to two-photon input. The analytical expression of the two- 
photon output state can be obtained if Eq. (17) can be calculated. Eq. 0 only 


contains the singular field operators and the Heisenberg-picture system operators which 
can be solved using their QSDEs. 

Interestingly, the four terms in Eq. 0 are not the most general form for the 


decomposition Eq. (16). We are able to further decompose {Li} as 


Li = eTL 


„ ^ , L * = l,2,---,i7, (18) 

with {6i} being constant column vectors. In other words, the coupling operator Li can 
be written as linear combination of a set of component operators {L'^, m = 1, 2, • • M}. 


Note that Li = is a special case for Eq. (18) but K and M may not be the same 


in general. Using Eq. (18), we can re-express Eq. ( [I^ as the linear combination of the 
following four terms 

Ol = {QQs\bj{T2)bi{Ti)b\ti)fi{ti)b\t2)i2{t2)\QQs), 

On = {00s\L'^{T2)bi{Ti)b'^ {ti)^i{ti)b'^ {t2)^2{t2)\00s), 

Oni = {00s\bj{T2)L'^{Tip{ti)^i{tip{t2)^2{t2)\00s), 

Oiv = {00s\L'^{T2)L'^{Tip{ti)^i{tip{t2)^2{t2)\00s) (19) 
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for n = 1, 2, • • - , M. In most cases, it is more convenient to deal with Eq. ( |I^ rather than 
Eq. (17). This is because the information of the linear combination (e.g. coefficients 


of the component operators) does not affect the analytical computability of the output 
state. As a consequence, component operators normally have simpler forms compared 
to {Tj}, and so applying QSDE analysis on component operators often leads to simpler 
conditions. As we will show in Section the coupling operator for a two-level atom 
could be L = yd?(T_,K > 0. However, the QSDE analysis of L' = cr_ is sufficient for 
proving the analytical computability of the output state. 

The response to a general two-photon input state can be analyzed using the same 
formalism, simply by replacing the uncorrelated state Eq. 0 with the general two- 
photon state in the derivations. For example, if the system couples to the environment 
via a single channel, a general two-photon input is written as 

/ oo noo 

dti / dt2i{ti,t2)h\ti)h\t2)\t)), ( 20 ) 

■OO J —oo 

with the normalization condition ff^dti dt 2 jf(ti,t 2 )l^dtidt 2 = 1. It is easy to see 


that the derivations remain the same if we use the general input state Eq. (20) instead 


3. Analytical computability of the output field state 


In this section we focus on the analytical computability of Eq. (19). It is easy to see 


that Oj can be exactly calculated using the commutation relation of the field operators 
{bi,bl}. The calculation of the rest terms in Eq. (19) relies on the Heisenberg-picture 
dynamics of the operator L'^{t), which is characterized by the QSDE 

i„(i) = C>{L'jt)) + lLHt),L'jmt) + b\t)lL'Jt),Lm, (21) 

where [L"'(t), L'^{t)] = [[L\{t), L'^{t)], ■ ■ ■,[L^j^{t), L'^{t)]] is a row vector of commutators. 


By Eq. (21), we can write the QSDE of L'(t) in a vector form 

L'{t) = C\L'{t)) + B{t)b{t) + b\t)C{t), (22) 

where £l(L'(t)) = [£l(L'^(t)) ■ ■ ■ C"'[L' j^[t))Y. The operator matrix B is defined by 
B = {[L^,L(„]}, i.e. the (m, j)-th entry of B is the commutator Similarly, the 

operator matrix C is defined hy C = {[L'^,Lj\}, i.e. the (m,j)-th entry of C is the 
commutator [L'^,Lj\. Eq. (22) implies 

(00,|L'(t) = (00,|(£t(L'(t)) + B{t)b{f)). (23) 

Therefore, if the following conditions 

£t(L'(()) = 2lL'((), (24) 

(00,|B = (00,|B, (25) 


hold with constant matrices A and H, then Eq. (23) is a solvable ordinary differential 


equation (ODE). By Eq. (24), each £l(L(„(t)) is required to be a linear combination of 
the component operators. Furthermore, it is required that A must be Hurwitz. Hurwitz 
means that the real parts of the eigenvalues of A are strictly negative. A Hurwitz matrix 
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A can remove the instantaneous response of the system and keep only the steady-state 
dynamics. In stability theory, A being Hurwitz is equivalent to the asymptotic stability 
of the linear system Eq. In most cases, since the energy is conserved in the overall 
system (see Assumption [^ , the photons will eventually leak to the fields and the system 
will be stabilized to its ground state. As a result, A being Hurwitz could be a natural 
property of the systems considered in this paper. 

According to Eq. (25), the elements of B = {bmj} are given by (00^1 [lJ, = 
(00s|&mp In other words, |00<j) is an eigenvector of the commutators {[{L'^y, Lj]}. 
Solving Eq. (23) and letting R —)■ —oo we have 


{00,|L'(i) = (00,1 / e-''‘-'^Bb(r)dr, 


(26) 


which is in the form of a convolution. Obviously, a frequency-domain relation naturally 
follows from the convolution. Here we have used the condition A being Hurwitz, so that 
the instantaneous term in the solution of the ODE converges to zero as to — cx). It 
is worth mentioning that the steady-state solution Eq. (26) only contains excitations in 
the field, which is consistent with our discussion in the last section. 


Alternatively, we can write Eq. (26) as 


(00,|L^(t) = (00, 


gm{t - r)b{r)dr, 


(27) 


where gm{t — r) is the m-th column of e^^^ Using Eq. (27), we can write On as 


On = ( 00 , 


gn{'r2 - r)6(r)dr6i(ri)6'^(ti)^i(ti)6'^(t2)6(^2)|00,), (28) 


which can be readily calculated using the commutation relations. Similarly, we can 
write Ojv as 


r-T2 


Oiv = 


gn{r2 -r)dr{ms\b{r)]y^{Ti)b\ti)fi{ti)b\t2)i2{t2)\fifis)- (29) 


By Eq. (29), Om and Oiv are analytically computable if 
(00,|6j(T)L(„(Tl)6'(^l)^l(^l)iI*(<2)^2(^2)|00•) 


(30) 


is analytically computable for arbitrary j and r. Now the question is whether we can 
express Eq. (30) using field operators only so that we can again apply the commutation 
relations. 


The formal integration of Eq. (22) leads to 

L (t) = [ e^^^~^\B{r)b{r) + b\r)C{r))dr. 


(31) 


into Eq. (ISOj), and again making use of the commutation relation. 


Substituting Eq. 

it is straightforward to show that the sufficient condition for Eq. (30) to be analytically 
computable is that the following terms 

(00,|&j(r)[(L^)'^(r),L|„(r)]6j(tg)|00,), (00,|[L^(r), L|„(r)]6j(ti)6](t2)|00,) (32) 

are analytically computable for arbitrary i,j,m,n. The analytical computability of 
the two terms is largely dependent on the commutation property of the component 
operators. We consider two typical cases: 
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In this case, Eq. (32) can be readily computed using the commutation relation. 
According to the above discussion, the output state is analytically computable. Systems 
that satisfy this condition will be discussed in Section F 


S.2. 1(4)'. ij is a nontrivial operator and =constant 


Since is a constant, we have (00s|[A^(r), L^(r)]6|(ti)6j(t2)|00s) = 0. When 

1(4)'. 41 is a nontrivial operator, we can perform integration on Eq. (|^ to obtain 

(00,|(),(T)|(4)'(>-),4('')]*i(<«)|00«) 

= (00.16, (t) r *e'"»l’-l{li'(s).l(4)'(i').4(s)l|i'(ii) 

J —oo 

+ 6'(s)[l(4)'(s),4(ii)].i(ii)l}i>I((,)|00.> 

= S{s - «,)(00,|6,(t)| f dse“"-(’-‘>[4(s), l(4)’(i'). 4(f)ll|00.> 

J —OO 

+ 5(s-r)(00,|[ / (ise'=™('-*)[[(L|^)^(s),L|„(s)],Lj(s)]6|(fg)|00,) (33) 

J —OO 

if the condition 


£'(1(4)'. 41) = <;™4(4)'. 41. s(c„„) < 0 (34) 

is satished. Now recall Eq. ( [2^ , where we have proven that the component operators 
sandwiched between single photon state and vacuum state are exactly calculable. As a 
result, Eq. (33) is analytically computable if 

[l(4)'.4l.ri = Ar'. 


(35) 


with A being a. K x M constant matrix. 

We have derived a set of algebraic conditions which are easily checkable. The 
component operators {L'^} in these conditions may not necessarily be the actual 
coupling operators between the system and the fields. Nevertheless, the actual coupling 
operators must be linear combinations of the component operators. Consequently, for 
the purpose of computing the two-photon output state, the key thing is to identify a set 
of component operators {L'^} which satisfy the sufficient conditions. In the subsequent 
sections, we will show that the exact form of the time-domain output states can be 
obtained for the systems whose modelling and time-domain calculation are difficult using 
conventional approaches. Particularly, we obtain the exact real-time two-photon output 
state for one-channel and two-channel scattering by a two-level emitter. The scattering 
of a two-level emitter is a problem of critical importance and has been extensively 
studied in the literature. 


4. Passive linear network 

We consider a linear network composed of single-mode harmonic oscillators which are 
coupled together via interconnection. This may refer to an optical network. An 
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W (b) 

Figure 3. (a) The cascade interconnection of Gi and G 2 - The cascaded system has 
single input and output channel, (b) There is a linear interaction between Gi and G 2 ■ 
The combined system has two input and two output channels. 


optical linear network can be used to process the information encoded in photons 
[29l [30] . Due to the application in linear quantum computation and photonic circuitry, 
the linear optical network has been extensively studied in the engineering community 
[231 El EH El]. It is well known that there are two basic types of interconnections for 
the construction of a network: cascade and direct interaction. Coherent feedback and 
other types of interconnections can be built up from the two basic types. Therefore, we 
just need to prove the analytical computability of the response of cascaded and directly 
coupled linear systems. 

A cascade connection of two single-mode open cavities with internal modes 01,02 
and resonant frequencies Ui,U 2 is depicted in Figure [^a). Suppose the coupling 
operators for the cavities are Li = ydfiOi and L 2 = 1 /^ 202 , ^ 1,^2 > 0, respectively. 
ki,K 2 are the decay rates of the cavities. The output signal of the hrst cavity is fed 
into the second cavity as the input. Employing the network theory, the cascaded 
system is described by the coupling operator L = ydfiOi -b ydf 202 and Hamiltonian 
Hq = ^o|oi -b ^0202 -b ^ ~ O 1 O 2 ) |32lElEim]. L is a linear combination of 
oi and 02 and so we can let L' = [oi 02 ]^. By Eq. Q we have the linear equations 


CHa2) 


+ 

—-(ia ;2 + ^2)^2 ~ 


^-a2, 


a/^1^2 

2 


Oi, 


(36) 


whose vector form is 


c\l') = al\ A = 


2 ^ + ^2) 


A is Hurwitz if the real part of A is negative dehnite. Since 

3?(A) = <0 



(37) 


( 38 ) 
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b(t) 










Figure 4. The optomechanical system has single input and output channel. The 
decay of the mechanical mode is negligible compared to the decay of the cavity mode. 


holds for any Ki,K 2 > 0, we can conclude that A is Hurwitz. Furthermore, since the 
commutators [aj,ai] = — 1 , [ai,ai] = 0 , [ 03 , 02 ] = — 1 , [ 02 , 02 ] = 0 are constants, the 
two-photon output state can be exactly calculated. 

The directly coupled system as plotted in Figure |^b) is described by the linear 
coupling operator L = [yd^oi ^/K2a2]'^ and Hamiltonian Hq = ^a\ai -|- ^ 0^02 -b 
7(0301 -b O 1 O 3 ) [201E]- 7 is the coupling strength. The two subsystems are coupled via 
a linear interaction term 7(0301 -b O 1 O 2 ). In this case, we can still let L' = [oi 02 ]^ as 
in the cascade case. We can obtain the linear equations 


C\ai) = --(iwi -b Ki)oi - i 702 , 
£1(02) = -^(iw2 + K-2)a2 - i^ai. 
The coefficient matrix of C^{L') = AL' is 


A 


— KicOi -b Ki) 
-i7 


-17 

-Uiu2 + K2) 


(39) 


(40) 


which is Hurwitz. Therefore, the output state is exactly calculable for the directly 
coupled system. 


5. Optomechanical system 

An optomechanical system may follow linearized dynamical equation under certain 
circumstances [OH [35]. The optomechanical system as shown in Figure]^ is composed 
of a linear cavity and a mechanical oscillator in interaction. The linearized system 
Hamiltonian is given by Hq = -b -b 7 (cl -b c)(al -b a) [Mj, where 7 can 

be made a real number, c is the cavity mode and a is the mechanical mode. The 
optomechanical system couples to the external held by the coupling operator L = ^/Kc. 
Letting L) = c, we can derive 

C\c) = -^{iuc +n)c - i'y{a + a^). (41) 

Although the coupling operator L contains c only, the generator of c is dependent on 
a, ah This observation motivates us to dehne L' = [c A a a^]^. Note that L is still 
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a linear combination of the component operators, with the coefficients on c\ a, being 
zero. Since 


£■^( 0 ) = --iuma - i7(c + c^), 
we can express CI{L') as a linear equation 


c\l') = AL' = 


( -\\ 


\ 


fuJc + k) 
0 

-iq 

17 


0 

-iq 

17 


-17 

-i7 

0 


-17 

-i7 

0 

I'Wm 


(42) 




/ 


(43) 


A is Hurwitz due to 3fJ(y4) < 0, which is also a consequence of the passivity of the system. 
Additionally, the commutations between c, c\ a, are all constants. Thus the response 
of the optomechanical system to two-photon input can be exactly computed. 


6. Two-level emitter 

We consider a two-level emitter in interaction with the photons propagating in the 
waveguide. The coupling of the two-level emitter to the optical helds is often modelled 
by where hiiuj) is the annihilation operator for the i-th 

mode of the held. When the photons propagate along the waveguide unidirectionally, 
there is only one coupling channel. If there are two modes for the travelling photons 
(e.g. left-propagating and right-propagating), we should model the interaction using two 
coupling channels. To be specihc, the system is coupled to the left-going mode via one 
channel, and coupled to the right-going mode via another channel. Since the interaction 
is energy-preserving, the system is passive when there exists no additional control that 
pumps energy into the system. We apply the results of Sec. |^to two processes of wide 
interests: stimulated emission of a two-level atom and the inelastic scattering of two 
photons. 

6.1. Simulation of stimulated emission 

The stimulated emission can be modelled by a two-level emitter interacting with a 
single-channel input held via L = yfK{a-. Without loss of generality we assume the 
system Hamiltonian is Hq = 0. In particular, introducing a free Hamiltonian of the 
form Hq = will only induce an additional harmonic component with frequency bJc 
in the pulse functions. 

When the system is fully excited and coupled to vacuum, the spontaneous emission 
rate is ki. However, if a second incoming photon interacts with the population-inverted 
emitter, the emission of a photon may either accelerate or slow down, depending on the 
exact form of the input pulse. 

It is easy to verify the commutation relations of L' = a- by [(j+,cr_] = 
cr^, [(T+, [(j+,cr_]] = —2(j+ and (00^1 [cr+, (j_] = (00s|(j^ = —(00s|. Also, we can obtain 

cT_(t) = - + \f^i(yz{t)b{t) 
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a^{t) = - + az{t))-2y/KS{h\t)a-{t) + a+{t)h{t)). (44) 

Therefore, we can conclude and (T_]) = C'^^az) = —ki(/ + cr^). 

Noting that the output field state is exactly solvable, we just need to design the 
pulse shapes so that the two-photon interaction could simulate the stimulated emission 
process. Intuitively, the first photon, which is followed immediately by the second 
photon, should be able to fully excite the system from the ground state to the excited 
state. For this reason, we choose the pulse function of the hrst photon to be the following 
form 


= -\/ 7 e “ ‘(1 


(45) 


where u{ti) is the Heaviside step function and 7 is a controllable parameter. Eq. (45) is 
the famous rising exponential pulse which can perfectly transfer the single photon to the 
two-level system at t = 0 [36l |37] when we let 7 = ki. As a result, the second incoming 
photon should be dehned on t 2 € (0,-|-cx)). Following the procedures in Sec. we can 
exactly calculate the output field state to be 




<iri<ir2^'(Ti.T2)5*(Ti)()t(r2)|00,), 


(46) 


with the analytical form of the pulse function given by 

+ K-i [ e~^^'"^~^'^^~''^f2{s)ds, T2 >ti> 0 . 


(47) 


It is straightforward to employ this analytical form to study the optimal stimulated 
emission, which could be used for photon amplification (TB]. For example, if the pulse 
function of the second photon is defined as ^ 2 (^ 2 ) = with H 2 being a 

controllable parameter, then we have 

e'(ri, T 2 ) = + (1 - (48) 

Ki — K,2 ~ f^2 

for r 2 > Ti > 0. Obviously the second photon cannot increase the decay rate if K 2 < 
and so we only consider the case K 2 > ki. We aim to maximize the component of 
Eq. (48) that decays with the rate K 2 , and minimize the component with the decay rate 
Ki- This is done by letting 

Kl — K2 

In this case, the faster decay rate K 2 dominates. The emission of the photon tends to 
synchronize with the faster decay rate of the second incoming photon, causing photon 
bunching in the output held. As a consequence, we will observe both photons in the 
output held earlier than the typical spontaneous emission time. The optimal K 2 is thus 
given by 


= 0 . 


(49) 


«2 = (50) 

Since ^'(ti, T 2 ) is a real function, we can conveniently plot ^'(ti, r 2 ) to compare the 
performance of diherent H 2 , see FigureWe let Ki = 1. The maximal photon bunching 
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4(t^,T2),K2=1.2 



5{t,,T2),K2=10 



1 2 


Figure 5. The value of ^ {ti,T 2 ) for different K 2 , calculated using Eq. (48). The 
result is symmetric with respect to ri = T 2 . ki = 1. K 2 = 3 (left, lower row) gives the 
optimal performance. When K 2 = 0.5 or K 2 = 1.2, the emitter decays slowly. When 
K 2 = 10, there is a high probability that the interval between the emission of the two 
photons is longer than 1. 


is observed when /t 2 = 3. When K 2 = 10, although there is higher probability that 
the stimulated emission would happen in a short time ('ri,r 2 < 0.5), the probability 
of observing a delayed second photon is also high compared to /t 2 = 3 case, e.g., 
(ti 0.1, r 2 = 1). In other words, there is still significant probability for detecting 
the photon anti-bunching. 

The previous works [T^ [TU] model the stimulated emission as a waveguide 
containing an incident photon in interaction with a two-level excited atom. The output 
states are then obtained using a real-space approach, that is, by solving the Schrodinger 
equation for two-photon wavefunctions. In contrast, our approach does not rely on a 
stationary-state expansion and so the calculation is more straightforward. In particular, 
the optimal stimulated emission has been studied in [11] using the two-time correlation 
function of the output. As we have mentioned before, the correlation analysis can be 
easily done since we can obtain the exact output state. As a matter of fact, if the decay 
constant of the atom is normalized to 1, then the optimal K 2 would be 3 according to 


the calculations in [T6|, which is consistent with our result Eq. (50). 
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In the above we have considered the optimal stimnlated emission. The spontaneous 


emission with decay rate Ki will be enhanced if Eq. (50) holds, and both photons will 
tend to decay with the faster rate K 2 > ki. Next we will discuss the suppression 
of spontaneous emission when H 2 < ki. Again considering Eq. (48), we still need 
to maximize the component that decays with the slower rate K 2 - However, since 


> K 2 > 0, Eq. (50) is not feasible and so the optimal suppression of spontaneous 


emission cannot be realized. In order to prolong the life time of the excited atom, we 
need to consider alternative pulse function for the second photon. For example, we 
could dehne the piecewise pulse function ^ 2 (^ 2 ) = \/>^ 2 /^ [0)^] ^md 
^ 2 (^ 2 ) = 0 elsewhere. Then we have 


/ A K1K2 

^(^ 1 , 7 - 2 )= _ e ^ 


^{(1 + 


2 a/^ 


y~ 


+ ( 1 - 


2 aAi ^^-^(r 2 -rl) 


-)e 


}(51) 


^eK 2 T _ X " K1-K2' ' K1-K2 

for T > r 2 > Ti > 0. The condition to minimize the spontaneous emission is given by 


1 - 


- ^2 


= 0 , 


or 


K2 = Ki - 2^/i^, 


(52) 


(63) 


which requires Ki > 4. When we choose K 2 to satisfy Eq. (53), the hrst emission at ti 
is induced by the reduced decay rate ki — K 2 , and the spontaneous emission with the 
decay rate ki is maximally suppressed before the emission of the second photon. Note 
that ^ 2 (^ 2 ) is similar to the rising exponential pulse function which is designed to excite 
the atom. 


6.2. Two-channel and one-channel scattering 

As we have mentioned, two coupling channels are used for the modelling of the left-going 
and right-going optical helds which are scattered at the emitter. The coupling operator 
is given by L = = [a/^i ^/K 2 ]'^c^_. Let L' = (T_ and we can exactly 

calculate the output held state. Suppose the two input photons are separated in two 
channels, with Q{t),i = 1,2 being the pulse functions of the input photons travelling in 
the i-th channel. The output held state is expressed as 

/ oo poo 

/ dTidT2[f'n{Ti,T2){b\y + ^'22{Ti,T2){blY + ^[2iu^^2)b\bl]\0). (54) 

■oo J —OO 

There are three diherent two-photon components in the output held, namely, the 
probability of two photons in the hrst channel, two photons in the second channel, 
and one in the hrst and one in the second. 

We consider the probability P of observing at least one photon in the hrst output 
channel 

/ OO poo 

/ dridr2(|^u(H,r2)P + |^i2(ri,r2)n. 

■OO J —OO 


( 55 ) 
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Figure 6. Left: The transmission probability of a two-level emitter. Right: The 
transmission probability of a single-mode linear cavity in response to the same two- 
photon input Eq. (57). The decay rate of the cavity mode is k, which is the same 
as the decay rate of the emitter. There exists a rapid transition of P for a two-level 
emitter when 71 < 0.5. 


For simplicity we assume Ki = K 2 = k which is the case for a waveguide system (the 
emitter couples to the optical fields with equal strength). The expressions for T 2 ) 

and ^ 12 ('Ti, T 2 ) with ti > T 2 are calculated as 

^n(n,T-2) =-K / -/t / 

J —00 J —00 

J — 00 J —CO 

f ' [ dr[ei(s) 6 (r) + 6 (r) 6 (s)], 

J —00 J —00 

^i2(n,T-2) =-K / e~'^^'^^~'^^i2{s)ii{Ti)ds - K j e~'^^'^^~"\i{s)i2{T2)ds 

J —CO J —CO 

+^2 f ' f ' +ei(s) 6 (r)] 

J — 00 J —CO 

/ e^^^^+^^ds f dr[^i(s)^2(r) + 6(^^)6(s)] + 6(n)6(T‘2)- (56) 

J —00 J —00 

Suppose the pulse functions of the two input photons are given by 

q(ii) = -VTre^'‘(l-«(ti)), fete) =-V^e^'yi - afe)), (57) 
with 7 i ,72 being controllable parameters. In this case, the relation between the 
transmission probability P and the parameters 71,72 is shown in Figure Also we have 
calculated the two-photon output state if the two-level emitter is replaced by a linear 
cavity. The first observation is that the responses of the linear and two-level system 
are quite distinct. The nonlinearity induced by the two-level emitter could significantly 
decrease the transmission probability in a two-photon process. This conversion of the 
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transmission behaviour is consistent with the previous hndings, e.g. in [7| IT^[T71 ITT] . 
For instance, for a single photon incident on the emitter, the transmission probability 
goes to zero for long pulses and goes towards 1 as the pulse width goes to 0. For 
the latter case, the photon is no longer on resonance with the emitter. However, the 
presence of a second photon drastically changes the above behavior. This is because the 
two-photon state has different resonant frequencies compared to single-photon state. A 
two-photon state is in resonance if its energy matches the two-excitation eigenstate of 
the system [3HIIH]. As a two-level emitter can store at most one quanta, the photons 
cannot enter the emitter simultaneously. This leads to the nonlinear behaviour which 
is different from the single photon case. In addition, since photons do not interact in 
linear systems, the switching property of a linear cavity is the same for single-photon 
and two-photon states. Secondly, when we change the values of 71 and 72 for a linear 
system, the corresponding variation of P is smooth. However for a two-level emitter, 
the variation of P exhibits nonlinear behaviour when 71 is small, which is a signature 
of strong photon-photon correlation ME]. To be specihc, if we £x a small 71 , the 
transmission probability P may fluctuate from the minimum value to the maximum 
value for a small variation of 72 . 

For the purpose of further illustrating the nonlinearity of two-photon interaction, 
next we consider the one-channel scattering by dehning L = y/Ka^. In the literature, 
the analytical study of the two-photon scattering often makes use of frequency-domain 
scattering analysis mm ITO] or diagrammatic approaches m- Following the previous 
derivations, the time-domain pulse function of the output state can be exactly 

calculated. Dehning the Fourier transform of f{t) as f{oj) = f(t)e~‘^'^‘^^^dt, the 

frequency-domain representation of ^'(Ti,r 2 ) is given by 


^ (a;i,a;2) = ^i(n;i)^2(w2) + 6('^2)6('^i) 

—l^\- -^-r + .. . „-r]['C2(i^l)'Cl('^2) + 'C2(i^2)'Cl(<^l)] 


I 27ra;ii | 27ra;2i 


27rca2i)(f + 27ri:n2i) 

4^2 r 

(f -F 27ra;2i)(f + 27ra;ii) 


7 T [^2 (1:^1) 6 (1^2) + ^2 (1:^2) 6(1^1)] 


+U2-t) + ^ 6 ( 1^1 + ^2 


■ I -h 27rri 


I 27rri' 


't)], 

(58) 


where 'Ci(')) 'C 2 (') are the pulse functions of the hrst and second photon, respectively. The 
last term in the above equation characterizes the inelastic scattering of two photons. To 
be more specihc, the output photons with frequencies oji and UJ 2 can be generated by 
a pair of incident photons with diherent frequencies r and ui + u:2 — t. This energy¬ 
preserving inelastic scattering property matches the previous hndings in mmin]- 
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7. Conclusion 

We have proposed a QSDE approach to model the two-photon scattering process for 
a general quantum system and calculate the time-domain response. We have studied 
only two cases in Section which enable the exact calculation of the two-photon output 
state. Nevertheless, there may exist quantum systems which allow exact analysis but do 
not belong to the two cases. For example, the response of a two-level system embedded 
in an open cavity or a linear cavity with Kerr nonlinearity has been analytically studied 
in lini |38]. The specihc system operators considered in these works do not satisfy 
the sufficient conditions proposed in Section However, it is easy to show that the 
coefficient terms of the output state can be exactly calculated based on a QSDE analysis. 

As a result, these systems are all amendable to the QSDE approach proposed in this 
paper. It will be interesting to investigate whether there exists a more general condition 
which implies the analytical computability of two-photon response for all these systems. 

As pointed out in [12] , a photon must exist as a pulse in both waveguide and free- 
space. The analytical approach proposed in this paper is thus directly applicable to 
temporal pulse shaping [22] , as compared to the previous frequency-domain approaches. 
Moreover, the QSDE approach is extremely powerful in modelling a network of quantum 
systems, which has also been demonstrated in Section Therefore, the analysis of 
two-photon response for a complicated quantum system could be benefited from this 
research. We expect that our results may be more general and directly applicable to the 
design of on-chip quantum circuit, in which the propagating photons could be scattered 
by linear and two-level components. For such practical systems, the QSDE approach 
and network theory may need certain extension in order to study arbitrary system 
parameters, relaxation and non-Markovian effects HU. 
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